################################################################################
# Replication code for "Knowing is Disturbing": Main Part        
# Date: 2025-09-29                                            
# Author: Danqi Guo & Genia Kostka                            
################################################################################

### clear history
rm(list = ls())

## set working directory to the folder where the data and code are stored
# setwd("your working directory")

### Load library, if not installed yet, please install them first using install.packages("package name")
library(tidyverse)
library(ltm)
library(coefplot)
library(ggpubr)
library(estimatr)
library(sjPlot)

# load raw data
data<-read.csv("data_cleaned_rep.csv",sep=",")

################################################################################
# Table 1. 3x3 factorial design (manually fitting the 3x3 matrix with the tabulation results beloew)
################################################################################

table(data$condition)
# or use the following code, 0 = no mention, 1 = personalized, 2 = public
table (data$Surveillance,data$Censorship)

################################################################################
#Figure 2. Treatment Effect of Digital Control on Emotions and Public Attitudes
################################################################################

## Generate DV1 emotional responses by conducting Principal Component Analysis (PCA)
df_emotion<- data[ , c('Sadness','Anger',"Fear","Disgust","DoNotCare","Helplessness",
                       "Secured","Happiness","Curiosity","Surprised")] 

# check out correlation matrix 
cor_emotion<-cor(df_emotion)
mean(cor(df_emotion))
# run PCA using base R function princomp()
pca_emotion<-princomp(df_emotion) 
# information output
names(pca_emotion)
summary(pca_emotion,loadings=TRUE) # provided accurate variance but no eigenvalue (only sd)
# create emotion variables based on PCA scores
pc_emotion = pca_emotion$scores
cor(pc_emotion) ## 5 components are independent from each other.
df_components_emotion<-data.frame(pc_emotion)
# correlation matrix of the components and the original variable
round(cor(df_emotion,pca_emotion$scores),3)
data$negative_emo<-pc_emotion[,"Comp.1"]
data$neutral_emo<-pc_emotion[,"Comp.3"]
data$positive_emo<-pc_emotion[,"Comp.2"]

## given the negative pca loading for positive and neutral emotoin, we reverse the pc of both.
positive_emo.df<-data.frame(data$positive_emo)
positive_emo.df<-positive_emo.df*(-1) # reverse positive emotion pca score by multiplying -1 
data$positive_emo_rv<-positive_emo.df$data.positive_emo
neutral_emo.df<-data.frame(data$neutral_emo)
neutral_emo.df<-neutral_emo.df*(-1) # reverse neutral emotion pca score by multiplying -1
data$neutral_emo_rv<-neutral_emo.df$data.neutral_emo

## Generate DV2 two different public attittudes
data$general_attitude_fic<-rowSums(data[,c("Attitude.I_1","Attitude.I_2","Attitude.I_3","Attitude.I_4")])
data$general_attitude_real<-rowSums(data[,c("Attitude.II_1","Attitude.II_2","Attitude.II_3","Attitude.II_4")])


#below using sub data set: 5 conditions, exluding the compound treatments 
# subsetting those groups that receive single treatment 
subset <- data[data$condition =="PersonalCensor" | data$condition == "control" |
                 data$condition=="PersonalSurv"|data$condition=="PublicCensor" |
                 data$condition =="PublicSurv", ]

mod_positive<-lm(positive_emo_rv~condition,data=subset)
mod_negative<-lm(negative_emo~condition,data=subset)
mod_neutral<-lm(neutral_emo_rv~condition,data=subset)
mod_attitude_fic<-lm(general_attitude_fic~condition,data=subset)
mod_attitude_real<-lm(general_attitude_real~condition,data=subset)

conditionlab<-c("Personal\nCensorship","Personal\nSurveillance",
                "Public\nCensorship", "Public\nSurveillance")

plot_effectemotion_subset<-multiplot(mod_positive,mod_neutral,mod_negative,intercept= FALSE,title="",dodgeHeight = 0.3,
                                     ylab = "Condition",xlab="Treatment Effect on Emotions", pointSize = 2.5,horizontal=TRUE)+theme_bw()+
  xlim(-0.8,0.8)+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+scale_y_discrete(labels= conditionlab)+
  scale_color_manual(name="emotion",values=c("darkred","grey","darkgreen"),labels = c("negative","neutral","positive"))+
  scale_shape_manual(name="emotion",labels = c("negative","neutral","positive"),values=c(16,17,18))


#political attitudes: ficitive and real life situation
plot_effectattitude_subset<-multiplot(mod_attitude_fic,mod_attitude_real,intercept= FALSE,title="",dodgeHeight = 0.2,
                                      ylab = "Condition",xlab="Treatment Effect on Political Attitudes", 
                                      pointSize = 2.5,horizontal = TRUE)+theme_bw()+
  xlim(-1.5,0.5)+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+scale_y_discrete(labels= conditionlab)+
  scale_color_manual(name="scenario",values=c("darksalmon","grey4"),labels = c("fictitious","real-life"))+
  scale_shape_manual(name="scenario",labels = c("fictitious","real-life"),values=c(22,25))


#below using sub data set: 4 compounded treatments 
# subsetting those groups that receive single treatment 
subset2 <- data[data$condition=="control"| data$condition =="Personalboth" | 
                  data$condition == "PersonalCensorPublicSurv" |
                  data$condition=="PersonalSurvPublicCensor"|data$condition=="PublicBoth" , ]

mod_positive<-lm(positive_emo_rv~condition,data=subset2)
mod_negative<-lm(negative_emo~condition,data=subset2)
mod_neutral<-lm(neutral_emo_rv~condition,data=subset2)
mod_attitude_fic<-lm(general_attitude_fic~condition,data=subset2)
mod_attitude_real<-lm(general_attitude_real~condition,data=subset2)

conditionlab2<-c("Personal\nBoth","Personal Cens.\nPublic Surv.",
                "Personal Surv.\nPublic Cens.", "Public\nBoth")

#emotions: positive and negative sentiment
plot_effectemotion_subset2<-multiplot(mod_positive,mod_neutral,mod_negative,intercept= FALSE,title="",dodgeHeight = 0.3,
                                      ylab = "Condition",xlab="", pointSize = 2.5,horizontal=TRUE)+theme_bw()+
  xlim(-0.8,0.8)+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+scale_y_discrete(labels= conditionlab2)+
  scale_color_manual(name="emotion",values=c("darkred","grey","darkgreen"),labels = c("negative","neutral","positive"))+
  scale_shape_manual(name="emotion",labels = c("negative","neutral","positive"),values=c(16,17,18))


#political attitudes: ficitive and real life situation
plot_effectattitude_subset2<-multiplot(mod_attitude_fic,mod_attitude_real,intercept= FALSE,title="",dodgeHeight = 0.2,
                                       ylab = "Condition",xlab="", 
                                       pointSize = 2.5,horizontal = TRUE)+theme_bw()+
  xlim(-1.5,0.5)+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+scale_y_discrete(labels= conditionlab2)+
  scale_color_manual(name="scenario",values=c("darksalmon","grey4"),labels = c("fictitious","real-life"))+
  scale_shape_manual(name="scenario",labels = c("fictitious","real-life"),values=c(22,25))

###combine two ate-plots together 
plot_emotion<-ggarrange(plot_effectemotion_subset, plot_effectemotion_subset2,
                        ncol = 2,  common.legend = TRUE, legend = "right", align = "hv")

plot_attitude<-ggarrange(plot_effectattitude_subset, plot_effectattitude_subset2,
                         ncol = 2,  common.legend = TRUE, legend = "right", align = "hv")

# Plotting Figure 2
fig_2<-ggarrange(plot_emotion, plot_attitude,
                     nrow = 2,  common.legend = FALSE, legend = "right", align = "hv")



################################################################################
# Figure 3 Mediation effect of emotions on attitudes 
################################################################################
# data preparation for Figure 3 
# First run mediation analysis 
# Covariates string
covariates_str <- paste(
  "Age", "Gender", "Education", "Ethnicity", "MonthlyIncome", "CCPMembership",
  "rural", "OverseaExperience_yes", "nationalism_nor", "liberalism_nor", "vpnsavvy",
  "conservation_nor", "openness_nor", "conscientiousness_nor", "neuroticism_nor",
  "extraversion_nor", "agreeableness_nor", "LifeSat_1",
  sep = " + "
)

# Build formulas OUTSIDE a function
form_med_positive_a <- as.formula(
  paste("positive_emo_rv ~ condition +", covariates_str)
)
form_med_positive_bc_fic <- as.formula(
  paste("general_attitude_fic ~ condition + positive_emo_rv +", covariates_str)
)
form_med_positive_bc_real <- as.formula(
  paste("general_attitude_real ~ condition + positive_emo_rv +", covariates_str)
)
form_med_negative_a <- as.formula(
  paste("negative_emo ~ condition +", covariates_str)
)
form_med_negative_bc_fic <- as.formula(
  paste("general_attitude_fic ~ condition + negative_emo +", covariates_str)
)
form_med_negative_bc_real <- as.formula(
  paste("general_attitude_real ~ condition + negative_emo +", covariates_str)
)

# Fit models directly
med_positive_a <- lm(form_med_positive_a, data = data)
med_positive_bc_fic <- lm(form_med_positive_bc_fic, data = data)
med_positive_bc_real <- lm(form_med_positive_bc_real, data = data)
med_negative_a <- lm(form_med_negative_a, data = data)
med_negative_bc_fic <- lm(form_med_negative_bc_fic, data = data)
med_negative_bc_real <- lm(form_med_negative_bc_real, data = data)

## store all 8 results of positive - real/fic 
unique(data$condition) 
treatment_arms <- setdiff(unique(data$condition), "control")
# Mediator models: M ~ condition + covariates
mediator_models <- list(
  positive_emo_rv = med_positive_a,
  negative_emo    = med_negative_a
)
# Outcome models: Y ~ condition + M + covariates
# Nested by outcome then mediator (so we can pick the correct Y-model for each M)
outcome_models <- list(
  fic  = list(
    positive_emo_rv = med_positive_bc_fic,
    negative_emo    = med_negative_bc_fic
  ),
  real = list(
    positive_emo_rv = med_positive_bc_real,
    negative_emo    = med_negative_bc_real
  )
)

# Function to run all combinations of mediation analyses
run_all_mediations <- function(med_models, out_models, arms,
                               treat_var = "condition",
                               control_value = "control",
                               boot = TRUE, sims = 1000) {
  
  if (!length(med_models)) stop("No mediator models supplied.")
  if (!length(out_models)) stop("No outcome models supplied.")
  if (!length(arms)) stop("No treatment arms supplied.")
  
  res <- list()
  
  for (med_name in names(med_models)) {
    model_m <- med_models[[med_name]]
    
    for (out_name in names(out_models)) {
      # outcome model must match this mediator
      model_y <- out_models[[out_name]][[med_name]]
      if (is.null(model_y)) {
        stop(sprintf("Missing outcome model for outcome='%s' & mediator='%s'.",
                     out_name, med_name))
      }
      
      for (arm in arms) {
        key <- paste("results", med_name, out_name, arm, sep = "_")
        res[[key]] <- mediation::mediate(
          model.m = model_m,
          model.y = model_y,
          treat   = treat_var,
          mediator = med_name,          # name of M variable in your data
          control.value = control_value,
          treat.value   = arm,
          boot = boot,
          sims = sims                   # note: sims (plural)
        )
      }
    }
  }
  res
}
# Run all mediations, takes around 15-20 minutates
results <- run_all_mediations(
  med_models = mediator_models,
  out_models = outcome_models,
  arms       = treatment_arms
)

# plotting the results of mediation analysis 
# Function to extract estimates and confidence intervals from mediate object
extract_results <- function(med_res, label) {
  data.frame(
    Estimate = c(med_res$d0, med_res$z0, med_res$tau.coef),
    LowerCI = c(med_res$d0.ci[1], med_res$z0.ci[1], med_res$tau.ci[1]),
    UpperCI = c(med_res$d0.ci[2], med_res$z0.ci[2], med_res$tau.ci[2]),
    Effect = factor(c("ACME", "ADE", "Total Effect"), levels = c("ACME", "ADE", "Total Effect")),
    Condition = label
  )
}

# Labels for each result set
labels <- paste("Result", 1:8)

# Extract results for each mediate object
results_list <-mapply(extract_results, results, labels, SIMPLIFY = FALSE)

# Combine all results into a single data frame
results_df <- bind_rows(results_list)
# change the value label of Label into treatment groups 
results_df <- results_df %>%
  mutate(Condition = case_when(
    Condition == "Result 1" ~ "Public Surv",
    Condition == "Result 2" ~ "Public Both",
    Condition == "Result 3" ~ "Personal. Surv & Public Censor",
    Condition == "Result 4" ~ "Personal. Censor & Public Surv",
    Condition == "Result 5" ~ "Personal. Both",
    Condition == "Result 6" ~ "Public Censor",
    Condition == "Result 7" ~ "Personal. Censor",
    Condition == "Result 8" ~ "Personal. Surv",
    TRUE ~ Condition  # Keeps the original value if no condition matches
  ))
# add new column to indicate mediator and outcome

# Create the new column values
mediator_outcome <- rep(c("positive_fic", "positive_real", "negative_fic", "negative_real"), each = 24)

# Add the new column
results_df$Category <- mediator_outcome

# Explicitly set the order of the factor levels for the Effect column
results_df$Effect <- factor(results_df$Effect, levels = c("Total Effect","ADE","ACME"))
levels(results_df$Effect)

# Plotting with dodge
position_dodge <- position_dodge(width = 0.5)

results_df$Category <- factor(results_df$Category,
                              levels = c("positive_fic",
                                         "positive_real",
                                         "negative_fic",
                                         "negative_real"))
# Figure 3
fig_3<-ggplot(results_df , aes(x = Estimate, y = Effect, color = Condition, shape = Condition)) +
  geom_point(size = 3, position = position_dodge, stroke = 0.5) +  # Increase point size and stroke
  geom_errorbar(aes(xmin = LowerCI, xmax = UpperCI), width = 0.5, position = position_dodge, size = 0.6) +  # Increase error bar width and size
  labs(x = "Estimate", y = "Effect", title = "", color = "Treatment", shape = "Treatment") +
  theme_bw() + 
  geom_vline(xintercept=0, linetype="dashed", 
             color = "black", size=0.5)+ # Facet by Category with free y-scales)
  theme(legend.title = element_text(size = 10)) +
  scale_color_brewer(palette = "Dark2") + # Use a different color palette for better distinction
  scale_shape_manual(values = 1:8) +
  theme(legend.position = "bottom") +   # remove background box) +
  facet_wrap(~ Category, ncol = 2,
             labeller = labeller(Category = c(
               "positive_fic" = "Positive emotions: attitudes (fictitious)",
               "positive_real"= "Positive emotions: attitudes (real-life)",
               "negative_fic" = "Negative emotions: attitudes (fictitious)",
               "negative_real"= "Negative emotions: attitudes (real-life)"
             )))


################################################################################
# Table 3. Main Effects of Digital Surveillance and Digital Censorship on Emotions and Public Attitudes
################################################################################

##convert two treatment dimensions (surveillance and censorship) into factor 

data$surveillance_fac <- factor(data$Surveillance, 
                                levels = c(0, 1, 2),
                                labels = c("No Mention",  "Personalized","Public"))
data$censorship_fac <- factor(data$Censorship, 
                              levels = c(0, 1, 2),
                              labels = c("No Mention", "Personalized", "Public"))

# Define outcomes
outcomes <- c("positive_emo_rv", "negative_emo", "neutral_emo_rv", 
              "general_attitude_fic", "general_attitude_real")

# Function to create models
fit_models <- function(outcome) {
  # Main effects formula
  formula_main <- as.formula(paste0(outcome, " ~ surveillance_fac + censorship_fac + ", covariates_str))
  # Interaction model
  formula_interact <- as.formula(paste0(outcome, " ~ surveillance_fac * censorship_fac + ", covariates_str))
  
  list(
    mod2 = lm_robust(formula_main, data = data), # no interaction
    mod4 = lm_robust(formula_interact, data = data) #with interaction
  ) 
}

# Fit all models
models <- lapply(outcomes, fit_models)
names(models) <- outcomes

# Flatten the nested list into a single list of models
all_models <- unlist(models, recursive = FALSE)

### models without control, with and without interaction Table D2 
tab_3<-tab_model(all_models$positive_emo_rv.mod2,all_models$positive_emo_rv.mod4,
                 all_models$negative_emo.mod2,all_models$negative_emo.mod4,
                 all_models$neutral_emo_rv.mod2,all_models$neutral_emo_rv.mod4,
                 all_models$general_attitude_fic.mod2,all_models$general_attitude_fic.mod4,
                 all_models$general_attitude_real.mod2,all_models$general_attitude_real.mod4,
                 dv.labels = c("Positive Emotion \n (main)","Positive Emotion \n (main+interactions)",
                                "Negative Emotion \n (main)","Negative Emotion \n (main+interactions)",
                                "Neutral Emotion \n (main)","Neutral Emotion b\n (main+interactions)",
                                "Attitude towards fictitious scenario \n (main)", "Attitude towards fictitious scenario \n (main+interactions)",
                                "Attitude towards real-life scenario \n (main)", "Attitude towards real-life scenario \n (main+interactions)"),
                 digits=3, p.style = "stars",show.ci=FALSE,show.se = TRUE,collapse.se = TRUE,
                 # <<<<<< Hide the controls
                 rm.terms = c(
                   "Age", "Gender", "Education", "Ethnicity", "MonthlyIncome", "CCPMembership",
                   "rural", "OverseaExperience_yes", "nationalism_nor", "liberalism_nor", "vpnsavvy",
                   "conservation_nor", "openness_nor", "conscientiousness_nor", "neuroticism_nor",
                   "extraversion_nor", "agreeableness_nor", "LifeSat_1"
                 ))

################################################################################
## ADDITIONAL: 
## I. calculate attitude difference in real life and hypothetical situation 
## II. calcualte the reliability of the attitude scale
################################################################################
# I. calculate attitude difference in real life and hypothetical situation
# calculate group-level means and difference A - B
group_means_att <- data %>%
  group_by(condition) %>%
  summarise(
    Mean_real = mean(general_attitude_real, na.rm = TRUE),
    Mean_hypo = mean(general_attitude_fic, na.rm = TRUE),
    Mean_Diff = mean(general_attitude_real - general_attitude_fic, na.rm = TRUE),
    SD_Diff = sd(general_attitude_real - general_attitude_fic, na.rm = TRUE),
    N = n(),
    .groups = 'drop'
  )

# paired t-test per group
t_test_results_att <- data %>%
  group_by(condition) %>%
  summarise(
    t_statistic = t.test(general_attitude_real, general_attitude_fic, paired = TRUE)$statistic,
    p_value = t.test(general_attitude_real, general_attitude_fic, paired = TRUE)$p.value,
    .groups = 'drop'
  )

# merge results for full summary
final_results_diff_att<- left_join(group_means_att, t_test_results_att, by = "condition")

# add bonferroni-adjuested p-values and SE
final_results_diff_att <- final_results_diff_att %>%
  mutate(
    SE = SD_Diff / sqrt(N),
    p_adj = p.adjust(p_value, method = "bonferroni"),
    Sig_Original = ifelse(p_value < 0.05, "Significant", "Not Significant"),
    Sig_Bonferroni = ifelse(p_adj < 0.05, "Significant", "Not Significant")
  )
# view the results
print(final_results_diff_att)
diff_att_df<-data.frame(final_results_diff_att)
###caculate the percentage increase in attitude from fictitious to real-life scenario
(sum(diff_att_df$Mean_Diff/diff_att_df$Mean_hypo)*100)/9

#II. calcualte the reliability of the attitude scale

attitude_df_fic <- data %>%
  dplyr::select(Attitude.I_1, Attitude.I_2, Attitude.I_3,Attitude.I_4)

attitude_df_real <- data %>%
  dplyr::select(Attitude.II_1, Attitude.II_2, Attitude.II_3,Attitude.II_4)
# Cronbach’s alpha
cronbach.alpha(attitude_df_fic,CI=TRUE)
cronbach.alpha(attitude_df_real,CI=TRUE)

################################################################################
#                       End of replication code 
################################################################################

